package com.radar.algorithm;

import cn.hutool.core.util.ArrayUtil;
import cn.hutool.log.Log;
import cn.hutool.log.LogFactory;
import com.radar.common.Common;
import com.radar.common.Const;
import com.radar.common.GisConstants;
import com.radar.common.ProductType;
import com.radar.exception.RadarCaclError;
import com.radar.io.RadarGridData;
import com.radar.utils.RefObj;
import org.meteoinfo.data.GridData;
import org.meteoinfo.global.Extent;
import org.meteoinfo.global.MIMath;
import org.meteoinfo.global.PointD;
import org.meteoinfo.projection.Reproject;
import org.meteoinfo.projection.info.ProjectionInfo;
import org.meteoinfo.projection.proj4j.CoordinateTransformFactory;
import ucar.unidata.geoloc.Earth;
import ucar.unidata.geoloc.projection.proj4.EquidistantAzimuthalProjection;
import wcontour.KDTree;

import java.util.*;
import java.util.stream.IntStream;

public class Cacl {

    static Log log = LogFactory.get();

    static final CoordinateTransformFactory ctFactory = new CoordinateTransformFactory();

    public static byte[][] gridppi(double[] azi, byte[][] z) {
        return gridppi(azi, z, ProductType.REF);
    }

    /**
     * Polar radar gridding
     * Less than four points use average, the rest use bilinear interpolation
     *
     * @param azi azimuth
     * @param z   polar radar
     * @return
     */
    public static byte[][] gridppi(double[] azi, byte[][] z, ProductType type) {
        RefObj<double[]> newazi = new RefObj<>(null);
        RefObj<byte[][]> newz = new RefObj<>(null);
        sortAzi(azi, z, newazi, newz);
        azi = newazi.getValue();
        z = newz.getValue();
        int step = 1;
        int bins = z[0].length;
        int Points = bins / step;
        byte[][] gridz = new byte[2 * Points][2 * Points];
        double[] finalAzi = azi;
        byte[][] finalZ = z;
        IntStream.range(0, 2 * Points)
                .parallel()
                .forEach(i -> {
                    for (int j = 0; j < 2 * Points; j++) {
                        gridz[i][j] = -128;
                        double bin0 = Math.sqrt(Math.abs(i - Points) * Math.abs(i - Points) + Math.abs(j - Points) * Math.abs(j - Points));
                        if (bin0 > bins - 1) {
                            continue;
                        }
                        double cz = Math.atan2((i - Points), (j - Points));
                        double azi0 = cz / (Math.PI / 180) + 90;
                        azi0 = azi0 > 360 ? azi0 - 360 : azi0;
                        azi0 = azi0 < 0 ? azi0 + 360 : azi0;
                        gridz[i][j] = MathCommon.Code(interpolatePolar(azi0, finalAzi, bin0, finalZ, type), type);
                    }
                });
        return gridz;
    }

    /**
     * Interpolate item polar to grid
     *
     * @param azi0 angle
     * @param azi  azis
     * @param bin0 distince
     * @param z    polar data
     * @return
     */
    public static double interpolatePolar(double azi0, double[] azi, double bin0, byte[][] z, ProductType type) {
        int bins = z[0].length;
        int index1 = 0;
        int records = azi.length;
        double zv = Double.MAX_VALUE;
        for (int m = 0; m < records; m++) {
            if (azi0 <= azi[m] && Math.abs(azi[m] - azi0) < zv) {
                zv = Math.abs(azi[m] - azi0);
                index1 = m;
            }
            if (azi[m] > azi0) {
                break;
            }
        }
        int index2 = index1 + 1;
        if (index1 != 0) {
            if (index1 == records - 1) {
                index2 = index1 - 1;
            } else {
                double a = Math.abs(azi[index1 - 1] - azi0);
                double b = Math.abs(azi[index1 + 1] - azi0);
                if (a < b) {
                    index2 = index1 - 1;
                }
            }
        } else {
            index2 = records - 1;
        }

        int bin1 = (int) Math.floor(bin0);
        int bin2 = bin1 + 1;
        if (bin0 == bins - 1) {
            bin2 = bin1 - 1;
        }
        // Less than four points use average, the rest use bilinear interpolation
        if (z[index1][bin1] == -128
                || z[index1][bin2] == -128
                || z[index2][bin1] == -128
                || z[index2][bin2] == -128) {
            double s = 0;
            int c = 0;
            if (z[index1][bin1] != -128) {
                s += MathCommon.Decode(z[index1][bin1], type);
                c++;
            }
            if (z[index1][bin2] != -128) {
                s += MathCommon.Decode(z[index1][bin2], type);
                c++;
            }
            if (z[index2][bin1] != -128) {
                s += MathCommon.Decode(z[index2][bin1], type);
                c++;
            }
            if (z[index2][bin2] != -128) {
                s += MathCommon.Decode(z[index2][bin2], type);
                c++;
            }
            if (c != 0) {
                return s / c;
            }
            return GisConstants.UNDEF;
        }
//        double z1 = z[index1][bin1] + (z[index1][bin2] - z[index1][bin1]) * (bin0 - bin1);
        double z1 = MathCommon.Decode(z[index1][bin1], type) + (MathCommon.Decode(z[index1][bin2], type) - MathCommon.Decode(z[index1][bin1], type)) * (bin0 - bin1);
//        double z2 = z[index2][bin1] + (z[index2][bin2] - z[index2][bin1]) * (bin0 - bin1);
        double z2 = MathCommon.Decode(z[index2][bin1], type) + (MathCommon.Decode(z[index2][bin2], type) - MathCommon.Decode(z[index2][bin1], type)) * (bin0 - bin1);
        double a1 = azi0 - azi[index1];
        double a2 = azi[index2] - azi[index1];
        if (a2 != 0) {
            return z1 + (z2 - z1) * a1 / a2;
        } else {
            return z1;
        }
    }

    /**
     * Transform longitude and latitude coordinates to gridded polar radar data
     * return lon lat project radar
     *
     * @param grid       grid radar
     * @param elevation  if 0 : not Calculated height
     * @param lon        center lon
     * @param lat        centos lat
     * @param resoultion resoultion km
     * @param step       lon lat default 0.01
     */
    public static RadarGridData project(byte[][] grid, double elevation, double lon, double lat, double resoultion, double step) {
        double radius = grid.length / 2 * resoultion;
        double xself = resoultion * step;
        double yself = resoultion * step;
        Extent extent = Common.getExtentRadius(lon, lat, grid.length / 2 * resoultion * 1000);
        double minLon = extent.minX;
        double maxLon = extent.maxX;
        double minLat = extent.minY;
        double maxLat = extent.maxY;
        int xnum = (int) ((maxLon - minLon) / xself + 1);
        int ynum = (int) ((maxLat - minLat) / xself + 1);
        byte[][] data = new byte[ynum][xnum];
        double[] xarray = new double[xnum];
        double[] yarray = new double[ynum];
        IntStream.range(0, ynum).forEach(i -> yarray[i] = minLat + i * yself);
        IntStream.range(0, xnum).forEach(i -> xarray[i] = minLon + i * xself);
        // "+ellps=WGS84 +proj=aeqd +lon_0=" + lon + " +lat_0=" + lat + " +x_0=0.0 +y_0=0.0 +units=km +no_defs"
        ProjectionInfo factory = ProjectionInfo.factory("+ellps=WGS84 +proj=aeqd +lon_0=" + lon + " +lat_0=" + lat + " +x_0=0.0 +y_0=0.0 +units=km +no_defs");
        EquidistantAzimuthalProjection equidistantAzimuthalProjection = new EquidistantAzimuthalProjection(lat, lon, 0, 0, Earth.DEFAULT);
//        MathCommon.fill(data);
        IntStream.range(0, ynum)
                .parallel()
                .forEach(i -> {
                    for (int j = 0; j < xnum; j++) {
//                            ProjectionPoint projectionPoint = equidistantAzimuthalProjection.latLonToProj(yarray[i], xarray[j]);
//                            PointD pointD = new PointD(projectionPoint.getX(), projectionPoint.getY());
                        PointD pointD = Reproject.reprojectPoint(xarray[j], yarray[i], GisConstants.lonlat, factory);
                        pointD.X = (int) ((pointD.X + radius) * Math.cos(elevation * Const.DEG2RAD) / resoultion);
                        pointD.Y = (int) ((pointD.Y + radius) * Math.cos(elevation * Const.DEG2RAD) / resoultion);
                        if (pointD.X < 0 || pointD.Y < 0 || pointD.X >= grid.length || pointD.Y >= grid.length) {
                            data[i][j] = -128;
                        } else {
                            data[i][j] = grid[(int) pointD.Y][(int) pointD.X];
                        }
                    }
                });
        RadarGridData gridData = new RadarGridData();
        gridData.xArray = xarray;
        gridData.yArray = yarray;
        gridData.missingValue = GisConstants.UNDEF;
        gridData.radarData = data;
        gridData.projInfo = GisConstants.lonlat;
        return gridData;
    }

    /**
     * Sort the data in order of azimuth from smallest to largest
     *
     * @param Azi
     * @param Z
     * @param newAzi
     * @param newZ
     */
    public static void sortAzi(double[] Azi, byte[][] Z, RefObj<double[]> newAzi, RefObj<byte[][]> newZ) {
        int N = 0;
        double a0 = Azi[Azi.length - 1];
        int rm = (int) Arrays.stream(Azi).filter(i -> MIMath.doubleEquals(GisConstants.UNDEF, i)).count();
        if (a0 != Arrays.stream(Azi).max().getAsDouble()) {
            for (int i = Azi.length - 2; i > 0; i--) {
                if (Azi[i] > a0 && Azi[i - 1] <= a0) {
                    N = i;
                    break;
                }
            }
        } else {
            double min = Arrays.stream(Azi).filter(i -> !MIMath.doubleEquals(GisConstants.UNDEF, i)).min().getAsDouble();
            N = ArrayUtil.indexOf(Azi, min);
        }
        if (N > Azi.length - N) {
            N = Azi.length - N;
        }
        double[] Azi2 = new double[Azi.length - N];
        byte[][] Z2 = new byte[Azi.length - N][Z[0].length];
        for (int i = 0; i < Azi.length - N; i++) {
            Azi2[i] = Azi[i + N];
            System.arraycopy(Z[i + N], 0, Z2[i], 0, Z[0].length);

        }

        double[] nAzi = Azi.clone();
        byte[][] nZ = Z.clone();
        for (int i = 0; i < Azi2.length; i++) {
            double min = Arrays.stream(Azi2).min().getAsDouble();
            int index = ArrayUtil.indexOf(Azi2, min);
            nAzi[i] = Azi2[index];
            for (int j = 0; j < Z2[0].length; j++) {
                nZ[i][j] = Z2[index][j];
            }
            Azi2[index] = Double.MAX_VALUE;
        }
        double[] rmAzi = new double[nAzi.length - rm];
        byte[][] rmnZ = new byte[nAzi.length - rm][nZ[0].length];
        for (int i = rm; i < nAzi.length; i++) {
            rmAzi[i - rm] = nAzi[i];
            rmnZ[i - rm] = nZ[i];
        }
        newAzi.setValue(rmAzi);
        newZ.setValue(rmnZ);
    }

    /**
     * Sort the data in descending order of elevation angle
     *
     * @param Ele
     * @param Z
     * @param newEle
     * @param newZ
     */
    public static void sortEle(double[] Ele, byte[][] Z, RefObj<double[]> newEle, RefObj<byte[][]> newZ) {
        int N1 = Arrays.binarySearch(Ele, Arrays.stream(Ele).min().getAsDouble());
        int N2 = Arrays.binarySearch(Ele, Arrays.stream(Ele).max().getAsDouble());
        if (N1 > N2) {
            int tmp = N2;
            N2 = N1;
            N1 = tmp;
        }
        double[] Azi2 = new double[N2 - N1 + 1];
        byte[][] Z2 = new byte[N2 - N1 + 1][Z[0].length];
        for (int i = 0; i < N2 - N1 + 1; i++) {
            Azi2[i] = Ele[i + N1];
            System.arraycopy(Z[i + N1], 0, Z2[i], 0, Z[0].length);
        }
        double[] nEle = Ele.clone();
        byte[][] nZ = Z.clone();
        for (int i = 0; i < Azi2.length; i++) {
            int index = ArrayUtil.indexOf(Azi2, Arrays.stream(Azi2).min().getAsDouble());
            nEle[i] = Azi2[index];
            for (int j = 0; j < Z2[0].length; j++) {
                nZ[i][j] = Z2[index][j];
            }
            Azi2[index] = Double.MAX_VALUE;
        }
        newEle.setValue(nEle);
        newZ.setValue(nZ);
    }

    /**
     * Composite Ref
     *
     * @param z
     * @return
     */
    public static byte[][] cr(byte[][][] z) {
        byte[][] cr = new byte[z[0].length][z[0][0].length];
        for (int y = 0; y < z[0].length; y++) {
            for (int x = 0; x < z[0][0].length; x++) {
                double min = GisConstants.UNDEF;
                byte b = -128;
                cr[y][x] = -128;
                for (int layer = 0; layer < z.length; layer++) {
                    if (z[layer][y][x] != -128) {
                        double c = MathCommon.DecodeRef(z[layer][y][x]);
                        if (c > min) {
                            min = c;
                            b = z[layer][y][x];
                        }
                    }
                }
                cr[y][x] = b;
            }
        }
        return cr;
    }

    public static byte[][] quickCR(byte[][][] Z, double[][] Azi, double[] Ele, double binLength, double maxZ) {
        // 库数
        int bins = Z[0][0].length;
        byte[][] cr = new byte[2 * bins + 1][2 * bins + 1];
        int x0 = bins;
        int y0 = bins;
        for (int i = 0; i < 2 * bins + 1; i++) {
            for (int j = 0; j < 2 * bins + 1; j++) {
                int x = (j - x0);
                int y = (i - y0);
                double a0 = Math.atan2(y, x) * 180 / Math.PI;
                if (a0 < 0) {
                    a0 += 360;
                }
                double v = GisConstants.UNDEF;
                byte c = -128;
                for (int m = 0; m < Z.length; m++) {
                    double l = Math.sqrt(x * x + y * y);
                    double bin = (l / Math.cos(Ele[m] * Math.PI / 180));
                    double z = Math.sqrt((Math.pow(bin, 2) - Math.pow(l, 2))) * binLength;
                    if (bin < Z[0][0].length && z < maxZ) {
                        // 计算方位角
                        int aziIndex = -1;
                        double azi = Double.MAX_VALUE;
                        for (int n = 0; n < Z[0].length; n++) {
                            if (Math.abs(Azi[m][n] - a0) < azi) {
                                azi = Math.abs(Azi[m][n] - a0);
                                aziIndex = n;
                            }
                        }
                        if (aziIndex != -1) {
                            if (v < MathCommon.DecodeRef(Z[m][aziIndex][(int) bin])) {
                                v = MathCommon.DecodeRef(Z[m][aziIndex][(int) bin]);
                                c = Z[m][aziIndex][(int) bin];
                            }
                        }
                    }
                }
                cr[i][j] = c;
            }
        }
        return cr;
    }

    /**
     * Echo Tops
     *
     * @param z          ppis z
     * @param et         base GridData
     * @param threshold  minimum value of reflectivity to be taken into calculation
     * @param elevation
     * @param radarHigh  height of radar
     * @param resoultion units km
     * @return
     */
    public static RadarGridData et(byte[][][] z, RadarGridData et, double threshold, double[] elevation, double radarHigh, double resoultion) {
        et.caclData = new short[et.getYNum()][et.getXNum()];
        double cenLon = et.xArray[et.getXNum() / 2];
        double cenLat = et.yArray[et.getYNum() / 2];
        for (int i = 0; i < et.getYNum(); i++) {
            for (int j = 0; j < et.getXNum(); j++) {
                et.caclData[i][j] = -128;
                for (int m = elevation.length - 1; m > 0; m--) {
                    double z0 = MathCommon.DecodeRef(z[m][i][j]);
                    if (z[m][i][j] != -128 && z0 > threshold) {
                        double distince = MathCommon.getDistanceKM(et.xArray[j], cenLon, et.yArray[i], cenLat);
                        et.caclData[i][j] = (short) ((Math.tan(elevation[m] * Const.DEG2RAD) * distince * resoultion + radarHigh / 1000.0) * 100.0);
                        break;
                    }
                }
            }
        }
        return et;
    }

    /**
     * Echo Bottom
     *
     * @param z          ppis z
     * @param eb         base GridData
     * @param threshold  minimum value of reflectivity to be taken into calculation
     * @param elevation
     * @param radarHigh  height of radar
     * @param resoultion units km
     * @return
     */
    public static RadarGridData eb(byte[][][] z, RadarGridData eb, double threshold, double[] elevation, double radarHigh, double resoultion) {
        eb.caclData = new short[eb.getYNum()][eb.getXNum()];
        double cenLon = eb.xArray[eb.getXNum() / 2];
        double cenLat = eb.yArray[eb.getYNum() / 2];
        for (int i = 0; i < eb.getYNum(); i++) {
            for (int j = 0; j < eb.getXNum(); j++) {
                eb.caclData[i][j] = -128;
                for (int m = 0; m < elevation.length; m++) {
                    double z0 = MathCommon.DecodeRef(z[m][i][j]);
                    if (z[m][i][j] != -128 && z0 > threshold) {
                        double distince = MathCommon.getDistanceKM(eb.xArray[j], cenLon, eb.yArray[i], cenLat);
                        eb.caclData[i][j] = (short) ((Math.tan(elevation[m] * Const.DEG2RAD) * distince * resoultion + radarHigh / 1000.0) * 100.0);
                        break;
                    }
                }
            }
        }
        return eb;
    }

    /**
     * Calculate vertically integrated liquid.
     *
     * @param z          ppis z
     * @param vil        base GridData
     * @param threshold  minimum value of reflectivity to be taken into calculation
     * @param elevation
     * @param resoultion units km
     * @return
     */
    public static RadarGridData vil(byte[][][] z, RadarGridData vil, double threshold, double[] elevation, double resoultion) {
        vil.caclData = new short[vil.getYNum()][vil.getXNum()];
        double vilConst = 3.44e-6;
        double cenLon = vil.xArray[vil.getXNum() / 2];
        double cenLat = vil.yArray[vil.getYNum() / 2];
        for (int i = 0; i < vil.getYNum(); i++) {
            for (int j = 0; j < vil.getXNum(); j++) {
                vil.caclData[i][j] = -128;
                double distince = MathCommon.getDistance(vil.xArray[j], cenLon, vil.yArray[i], cenLat) * resoultion;
                double hi = Math.tan(0.99 / 2 * Const.DEG2RAD) * distince;
                double m1 = 0;
                int minIndex = -1;
                int maxIndex = -1;
                int m;
                for (m = 0; m < elevation.length - 1; m++) {
                    double ht = distince * (Math.sin(elevation[m + 1] * Const.DEG2RAD) - Math.sin(elevation[m] * Const.DEG2RAD));
                    if (MathCommon.DecodeRef(z[m][i][j]) > threshold) {
                        if (minIndex == -1) {
                            minIndex = m;
                        }
                        if (maxIndex == -1) {
                            maxIndex = m;
                        }
                    }
                    if (MathCommon.DecodeRef(z[m + 1][i][j]) > threshold) {
                        maxIndex = m + 1;
                    }
                    if (MathCommon.DecodeRef(z[m][i][j]) > threshold && MathCommon.DecodeRef(z[m + 1][i][j]) > threshold) {
                        double factor = Math.pow(((r2z(MathCommon.DecodeRef(z[m][i][j])) + r2z(MathCommon.DecodeRef(z[m + 1][i][j]))) / 2), (4.0 / 7.0));
                        m1 += vilConst * factor * ht;
                    }
                }
                if (minIndex == -1) {
                    continue;
                }
                double mb = vilConst * Math.pow(r2z(MathCommon.DecodeRef(z[minIndex][i][j])), (4.0 / 7.0)) * hi;
                double mt = vilConst * Math.pow(r2z(MathCommon.DecodeRef(z[m][i][j])), (4.0 / 7.0)) * hi;
                vil.caclData[i][j] = (short) ((m1 + mb + mt) * 100);
            }
        }
        return vil;
    }

    public static double r2z(double z) {
        return Math.pow(10, z / 10.0);
    }

    /**
     * cappi Single layer height bilinear interpolation
     *
     * @param azis
     * @param ele
     * @param alt
     * @param z
     * @return
     */
    public static byte[][] heightInterpolation(List<double[]> azis, double[] ele, double alt, List<byte[][]> z, double resoultion, ProductType type) {
        alt = alt / resoultion;
        int bins = z.get(0)[0].length;
        byte[][] data = new byte[bins * 2][bins * 2];
        for (int i = 0; i < data.length; i++) {
            for (int j = 0; j < data[0].length; j++) {
                double angle = MathCommon.getAngleXY(i, j, data.length / 2, data.length / 2);
                double distince = Math.sqrt(Math.pow(i - data.length / 2, 2) + Math.pow(j - data[0].length / 2, 2));
                data[data.length - 1 - i][j] = MathCommon.Code(bilinearInterpolat(angle, azis, ele, alt, distince, z, type), type);
            }
        }
        return data;
    }

    /**
     * cappi Single layer height bilinear interpolation by lon lat
     *
     * @param azis
     * @param ele
     * @param alt
     * @param z
     * @param resoultion
     * @param type
     * @param xarray
     * @param yarray
     * @param centerLat
     * @param centerLon
     * @return
     */
    public static byte[][] heightInterpolation(List<double[]> azis, double[] ele, double alt, List<byte[][]> z, double resoultion, ProductType type, double[] xarray, double[] yarray, double centerLat, double centerLon) {
        alt = alt / resoultion;
        byte[][] data = new byte[yarray.length][xarray.length];
        for (int i = 0; i < data.length; i++) {
            for (int j = 0; j < data[0].length; j++) {
                double angle = MathCommon.getAngleXY(yarray[i], xarray[j], centerLat, centerLon);
                double distince = MathCommon.getDistance(xarray[j], centerLon, yarray[i], centerLat) / resoultion;
                data[data.length - 1 - i][j] = MathCommon.Code(bilinearInterpolat(angle, azis, ele, alt, distince, z, type), type);
            }
        }
        return data;
    }


    public static double bilinearInterpolat(double angle, List<double[]> azis, double[] ele, double alt, double distince, List<byte[][]> z, ProductType productType) {
        double zv = Double.MAX_VALUE;
        // find ele index
        double d = Math.sqrt(Math.pow(distince, 2) + Math.pow(alt, 2));
        double e1 = 90 - Math.atan2(d, alt) * 180 / Math.PI;
        int bin = (int) Math.floor(distince);
        int ele_index_0 = -1;
        int ele_index_1 = -1;
        for (int i = 0; i < ele.length - 1; i++) {
            if (e1 >= ele[i] && e1 < ele[i + 1]) {
                ele_index_0 = i;
                ele_index_1 = i + 1;
                break;
            }
        }
        if (ele_index_0 == -1) {
            return GisConstants.UNDEF;
        }
        int ele0_angle_0 = -1;
        int ele1_angle_0 = -1;
        // find ele0 angle
        double[] azi = azis.get(ele_index_0);
        int records = azi.length;
        for (int m = 0; m < records; m++) {
            if (angle <= azi[m] && Math.abs(azi[m] - angle) < zv) {
                zv = Math.abs(azi[m] - angle);
                ele0_angle_0 = m;
            }
            if (azi[m] > angle) {
                break;
            }
        }
        if (zv > 2) {
            return GisConstants.UNDEF;
        }
        zv = Double.MAX_VALUE;
        int ele0_angle_1 = ele0_angle_0 + 1;
        if (ele0_angle_0 != 0) {
            if (ele0_angle_0 == records - 1) {
                ele0_angle_1 = ele0_angle_0 - 1;
            } else {
                double a = Math.abs(azi[ele0_angle_0 - 1] - angle);
                double b = Math.abs(azi[ele0_angle_0 + 1] - angle);
                if (a < b) {
                    ele0_angle_1 = ele0_angle_0 - 1;
                }
            }
        } else {
            ele0_angle_1 = records - 1;
        }
        // find ele1 angle
        azi = azis.get(ele_index_1);
        records = azi.length;

        for (int m = 0; m < records; m++) {
            if (angle <= azi[m] && Math.abs(azi[m] - angle) < zv) {
                zv = Math.abs(azi[m] - angle);
                ele1_angle_0 = m;
            }
            if (azi[m] > angle) {
                break;
            }
        }
        if (zv > 2) {
            return GisConstants.UNDEF;
        }
        int ele1_angle_1 = ele1_angle_0 + 1;
        if (ele1_angle_0 != 0) {
            if (ele1_angle_0 == records - 1) {
                ele1_angle_1 = ele1_angle_0 - 1;
            } else {
                double a = Math.abs(azi[ele1_angle_0 - 1] - angle);
                double b = Math.abs(azi[ele1_angle_0 + 1] - angle);
                if (a < b) {
                    ele1_angle_1 = ele1_angle_0 - 1;
                }
            }
        } else {
            ele1_angle_1 = records - 1;
        }
        try {
            if (bin >= z.get(0)[0].length) {
                return GisConstants.UNDEF;
            }
            double v1 = MathCommon.Decode(z.get(ele_index_0)[ele0_angle_0][bin], productType);
            double v2 = MathCommon.Decode(z.get(ele_index_0)[ele0_angle_1][bin], productType);
            double v3 = MathCommon.Decode(z.get(ele_index_1)[ele1_angle_0][bin], productType);
            double v4 = MathCommon.Decode(z.get(ele_index_1)[ele1_angle_1][bin], productType);
            if (MIMath.doubleEquals(v1, GisConstants.UNDEF) || MIMath.doubleEquals(v2, GisConstants.UNDEF) || MIMath.doubleEquals(v3, GisConstants.UNDEF) || MIMath.doubleEquals(v4, GisConstants.UNDEF)) {
                return GisConstants.UNDEF;
            }
            double linev1 = MathCommon.linearInterpolate(v1, v2, azis.get(ele_index_0)[ele0_angle_0], azis.get(ele_index_0)[ele0_angle_1], angle);
            double linev2 = MathCommon.linearInterpolate(v3, v4, azis.get(ele_index_1)[ele1_angle_0], azis.get(ele_index_1)[ele1_angle_1], angle);
            return MathCommon.linearInterpolate(linev1, linev2, ele[ele_index_0], ele[ele_index_1], e1);
        } catch (Exception e) {
            return GisConstants.UNDEF;
        }
    }

    /**
     * 孤立点消噪质控
     *
     * @param var
     * @return
     */
    public static double[][] nineMedianFilter(double[][] var) {
        return nineMedianFilter(var, 3);
    }

    /**
     * 孤立点消噪质控
     *
     * @param var
     * @return
     */
    public static double[][] nineMedianFilter(double[][] var, int z) {
        int rows = var.length;
        int columns = var[0].length;
        double[][] nvar = new double[rows][columns];
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < columns; j++) {
                nvar[i][j] = var[i][j];
            }
        }
        int zindex = z - 1;
        double 降噪阈值 = 0.75;
        double[] value = new double[z * z];
        for (int i = 1; i < rows - zindex; i++) {
            for (int j = 1; j < columns - zindex; j++) {
                int index = 0;
                for (int m = i - 1; m < i + zindex; m++) {
                    for (int n = j - 1; n < j + zindex; n++) {
                        value[index] = var[m][n];
                        index++;
                    }
                }
                List<Double> vv = new ArrayList<>();
                Arrays.stream(value).forEach(c -> vv.add(c));
                int cc = vv.size();
                for (int ii = 0; ii < cc; ii++) {
                    if (vv.get(ii) == GisConstants.UNDEF) {
                        vv.remove(ii);
                        ii--;
                        cc--;
                    }
                }
                if (vv.size() < (z * z * 降噪阈值)) {
                    nvar[i][j] = GisConstants.UNDEF;
                } else {
                    DoubleSummaryStatistics doubleSummaryStatistics = vv.stream().mapToDouble(Double::doubleValue).summaryStatistics();
                    nvar[i][j] = doubleSummaryStatistics.getAverage();
                }
            }
        }

        return nvar;
    }

    //region Deprecated keep code

    /**
     * project ccrs AzimuthalEquidistant to lon lat
     *
     * @param z
     * @param lons
     * @param lats
     * @param resoultion units km
     * @return
     * @throws RadarCaclError
     */
    @Deprecated
    public static GridData projectData(double[][] z, double[][] lons, double[][] lats, double resoultion, double lon, double lat, double elevation) {
        double xself = resoultion * 0.01;
        double yself = resoultion * 0.01;
        // max min value
        Extent extent = Common.getExtentRadius(lon, lat, z.length / 2 * resoultion * 1000);
        double minLon = extent.minX;
        double maxLon = extent.maxX;
        double minLat = extent.minY;
        double maxLat = extent.maxY;
        int xnum = (int) ((maxLon - minLon) / xself + 1);
        int ynum = (int) ((maxLat - minLat) / xself + 1);
        double[][] data = new double[ynum][xnum];
        double[] xarray = new double[xnum];
        double[] yarray = new double[ynum];
        IntStream.range(0, ynum).forEach(i -> yarray[i] = minLat + i * yself);
        IntStream.range(0, xnum).forEach(i -> xarray[i] = minLon + i * xself);
        // Because the speed is too slow, the kdTree reprojection is abandoned
        // projectKDTree(lons,lats,z,data,xself,xarray,yarray);
        // Because idw interpolation loses too much precision Use nearest point interpolation
        // project(z, data, xarray, yarray, lon, lat, resoultion);
        for (int i = 0; i < yarray.length; i++) {
            for (int j = 0; j < xarray.length; j++) {
                data[i][j] = GisConstants.UNDEF;
                double findx = xarray[j];
                double findy = yarray[i];
//                double c = Math.atan2((findy - lat), (findx - lon));
//                double angle = 360 - c / Const.DEG2RAD + 90;
//                angle = angle > 360 ? angle - 360 : angle;
                double angle = MathCommon.getAngle(findx, lon, findy, lat);
                double distance = MathCommon.getDistance(findx, lon, findy, lat);
                double deltav = Math.cos(angle * Const.DEG2RAD) * distance / Math.cos(elevation * Const.DEG2RAD);
                double deltah = Math.sin(angle * Const.DEG2RAD) * distance / Math.cos(elevation * Const.DEG2RAD);
                int yindex = (int) Math.floor(deltav / (resoultion * 1000.0) + z.length / 2);
                int xindex = (int) Math.floor(deltah / (resoultion * 1000.0) + z.length / 2);
                if (xindex > 0 && yindex > 0 && xindex < z[0].length && yindex < z.length && !MIMath.doubleEquals(z[yindex][xindex], GisConstants.UNDEF)) {
                    data[i][j] = z[yindex][xindex];
                }
            }
        }
        GridData gridData = new GridData();
        gridData.xArray = xarray;
        gridData.yArray = yarray;
        gridData.data = data;
        gridData.projInfo = GisConstants.lonlat;
        gridData.missingValue = GisConstants.UNDEF;
        return gridData;
    }

    /**
     * project lonlat by Azimuth idw
     * idw interpolation loses too much precision Use nearest point interpolation
     *
     * @param z
     * @param data
     * @param xarray
     * @param yarray
     * @param lon
     * @param lat
     * @param resoultion
     */
    @Deprecated
    public static void project(double[][] z, double[][] data, double[] xarray, double[] yarray, double lon, double lat, double resoultion, double elevation) {
        for (int i = 0; i < yarray.length; i++) {
            for (int j = 0; j < xarray.length; j++) {
                data[i][j] = GisConstants.UNDEF;
                double findx = xarray[j];
                double findy = yarray[i];
                double c = Math.atan2((findy - lat), (findx - lon));
                double angle = 360 - c / Const.DEG2RAD;
                angle = angle > 360 ? angle - 360 : angle;
                double distance = MathCommon.getDistance(findx, lon, findy, lat);
                double deltav = Math.cos(angle * Const.DEG2RAD) * distance * Math.cos(elevation * Const.DEG2RAD);
                double deltah = Math.sin(angle * Const.DEG2RAD) * distance * Math.cos(elevation * Const.DEG2RAD);
                int yindex = (int) Math.floor(deltav / (resoultion * 1000.0) + z.length / 2);
                int xindex = (int) Math.floor(deltah / (resoultion * 1000.0) + z.length / 2);
                if (xindex > 0 && yindex > 0 && xindex < z[0].length && yindex < z.length && !MIMath.doubleEquals(z[yindex][xindex], GisConstants.UNDEF)) {
                    data[i][j] = z[yindex][xindex];
                }
                // idw interpolation loses too much precision Use nearest point interpolation
//                double w = 0;
//                double mw = 0;
//                double value = 0;
//                if (xindex > 0 && yindex > 0 && xindex < z[0].length && yindex < z.length && !MIMath.doubleEquals(z[yindex][xindex], GisConstants.UNDEF)) {
//                    w = 1.0 / MathCommon.getDistanceKM(lons[yindex][xindex], findx, lats[yindex][xindex], findy);
//                    mw += w;
//                    value += z[yindex][xindex] * w;
//                }
//                if (xindex > 0 && yindex + 1 > 0 && xindex < z[0].length && yindex + 1 < z.length && !MIMath.doubleEquals(z[yindex + 1][xindex], GisConstants.UNDEF)) {
//                    w = 1.0 / MathCommon.getDistanceKM(lons[yindex + 1][xindex], findx, lats[yindex + 1][xindex], findy);
//                    mw += w;
//                    value += z[yindex + 1][xindex] * w;
//                }
//                if (xindex + 1 > 0 && yindex > 0 && xindex + 1 < z[0].length && yindex < z.length && !MIMath.doubleEquals(z[yindex][xindex + 1], GisConstants.UNDEF)) {
//                    w = 1.0 / MathCommon.getDistanceKM(lons[yindex][xindex + 1], findx, lats[yindex][xindex + 1], findy);
//                    mw += w;
//                    value += z[yindex][xindex + 1] * w;
//                }
//                if (xindex + 1 > 0 && yindex + 1 > 0 && xindex + 1 < z[0].length && yindex + 1 < z.length && !MIMath.doubleEquals(z[yindex + 1][xindex + 1], GisConstants.UNDEF)) {
//                    w = 1.0 / MathCommon.getDistanceKM(lons[yindex + 1][xindex + 1], findx, lats[yindex + 1][xindex + 1], findy);
//                    mw += w;
//                    value += z[yindex + 1][xindex + 1] * w;
//                }
//                if (w != 0) {
//                    data[i][j] = value / mw;
//                }
            }
        }
    }

    /**
     * KDTree project
     *
     * @param lons       source lons
     * @param lats       source lats
     * @param z          source z
     * @param data       target data
     * @param resoultion Search radius
     * @param xarray     target xarray
     * @param yarray     target yarray
     */
    public static void projectKDTree(double[][] lons, double[][] lats, double[][] z, double[][] data, double resoultion, double[] xarray, double[] yarray) {
        KDTree.Euclidean<double[]> kdTree = new KDTree.Euclidean(2);
        for (int i = 0; i < lats.length; i++) {
            for (int j = 0; j < lons.length; j++) {
                kdTree.addPoint(new double[]{lons[i][j], lats[i][j]}, new double[]{lons[i][j], lats[i][j], z[i][j]});
            }
        }
        for (int i = 0; i < data.length; i++) {
            for (int j = 0; j < data[0].length; j++) {
                double self = resoultion;
                int step = 0;
                data[i][j] = GisConstants.UNDEF;
                // Iterate 3
                while (step < 3) {
                    double findx = xarray[j];
                    double findy = yarray[i];
                    ArrayList<double[]> doubles = kdTree.ballSearch(new double[]{findx, findy}, self);
                    if (doubles.size() > 4) {
                        doubles.sort(Comparator.comparingDouble(v -> MathCommon.getDistance(v[0], findx, v[1], findy)));
                        double w = 0;
                        double value = 0;
                        for (int x = 0; x < 4; x++) {
                            if (!MIMath.doubleEquals(doubles.get(x)[2], GisConstants.UNDEF)) {
                                w += 1.0 / MathCommon.getDistance(doubles.get(x)[0], findx, doubles.get(x)[1], findy);
                                value += doubles.get(x)[2];
                            }
                        }
                        if (w != 0) {
                            data[i][j] = value / w;
                        }
                        break;
                    }
                    self *= 2;
                    step++;
                }
            }
        }
    }

    //endregion
}
